The maximum mass and radius of neutron stars and the nuclear symmetry energy 
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We calculate the equation of state of neutron matter with realistic two- and three-nucleon inter- 
actions using quantum Monte Carlo techniques, and illustrate that the short-range three-neutron 
interaction determines the correlation between neutron matter energy at nuclear saturation density 
and higher densities relevant to neutron stars. Our model also makes an experimentally testable 
prediction for the correlation between the nuclear symmetry energy and its density dependence - 
determined solely by the strength of the short-range terms in the three neutron force. The same 
force provides a significant constraint on the maximum mass and radius of neutron stars. 
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Since their discovery, neutron stars have remained our 
sole laboratory to study matter at supra-nuclear density 
and relatively low temperature. The equation of state 
(EoS) of matter at these densities is largely unknown 
but uniquely determines the structure of neutron stars 
and the relation between their mass (M) and radius (i?) . 
Matter that can support large pressure for a given en- 
ergy density (typically called a stiff EoS) will favor large 
neutron star radii for a given mass. Such an EoS also 
predicts large values for the maximum mass of a neutron 
star that is stable with respect to gravitational collapse 
to a black hole. Conversely, a high density phase that 
predicts a smaller pressure will result in more compact 
neutron stars and smaller maximum masses. 

The recent accurate measurement of a large neutron 
star mass M = 1.97 ± 0.04Msoiar in the system J1614- 
2230 provides strong evidence that the high density equa- 
tion of state is stiff [T]. Interestingly, attempts to infer 
neutron star radii have favored relatively small values 
ranging from 9 to 12 km [5H1]. Although the radius 
inference depends on specific model assumptions, these 
smaller radii imply a soft EoS in the vicinity of nuclear 
saturation density. Taken together, they indicate that 
the EoS of dense matter makes a transition from soft to 
stiff at supra-nuclear density. In this Rapid Communi- 
cation we show that the three-neutron force (3n) is the 
key microscopic ingredient that determines the nature of 
this transition. 

The importance of three-body forces in nuclear physics 
is well known, and quantum Monte Carlo (QMC) cal- 
culations of light nuclei have clarified its structure and 
strength. However, in these systems the dominant three- 
body force acts between two neutrons and proton or be- 
tween two protons and a neutron. While the force among 
three neutrons is important in light neutron-rich nuclei, 
the short distance behavior is not easily accessible [5]. 
Properties of large neutron-rich nuclei are potentially 
sensitive to this interaction, especially if the symmetry 
energy provides a reliable measure of the energy differ- 
ence between pure neutron matter and symmetric nu- 
clear matter at saturation density. There has been much 
recent progress in both theory and experiments to mea- 
sure the symmetry energy and its density dependence. 



as reviewed in Refs. [51 [7] ■ The symmetry energy is ex- 
pected to be in the range 32 ± 2 MeV. We explore this 
experimentally suggested range for the nuclear symme- 
try energy and show that a more precise determination 
is needed to adequately constrain the 3n interaction. 

In this work we solve the non-perturbative many-body 
nuclear Hamiltonian using the auxiliary field diffusion 
Monte Carlo (AFDMC) [8] method. Its accuracy in 
studying nuclear systems has been tested in light nu- 
clei [9j. The extension to include three-body forces 
in pure neutron rich systems is straightforward with 
no additional approximations within the AFDMC tech- 
nique |10j . and a comparison with the Green's function 
Monte Carlo (GFMC) has been extensively tested in neu- 
tron drops |TT] . We present results for the EoS of neutron 
matter using phenomenological two-neutron (2n) poten- 
tials, which provide an accurate description of nucleon- 
nucleon scattering data up to high energies, and study 
the role of the poorly constrained 3n interaction. 

In earlier work it has been established that the EoS 
in the density regime (1 — 3)po plays an essential role in 
determining the neutron star radius |12| . In this density 
regime, the 3n interaction plays a critical role because 
of a large cancellation between the attractive and repul- 
sive parts of the 2n interaction arising from the long and 
short distance behavior, respectively. Consequently, we 
find that the neutron star radius for a canonical mass 
of 1.4 Msoiar is especially sensitive to the 3n interaction. 
Although matter in the neutron star will contain a small 
admixture of protons, here we calculate the EoS of pure 
neutron matter for the following reasons. First, the struc- 
ture of the interactions between neutrons is simpler than 
those between neutron and protons. Second, these sim- 
pler interactions are amenable to QMC methods to solve 
the many-body problem as it is devoid of the complexities 
of the isospin dependent spin-orbit and three-nucleon po- 
tentials, and clustering effects likely in systems with pro- 
tons. Third, the fraction of protons required to ensure 
stability is small and is typically less than 10%. Finally, 
since generically neutron matter has higher pressure than 
matter containing any fraction of protons or strangeness 
in the form of hyperons or kaons, our results provide 
stringent upper bounds on the neutron maximum mass 
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and radius. 

To compute the EoS for neutron stars it is necessary 
to describe the nucleon-nucleon interactions at short dis- 
tances or large relative momenta up to p ~ '2pFn — 
660 MeV(p/po)^^^i where ppn is the Fermi momentum, 
p is the typical density in the neutron star core, and 
Po = 0.16 fm~'^ is the nuclear saturation density. Rela- 
tive momenta up to ppn are required in even a mean-field 
(Fermi gas) description, and the nn interaction scatters 
nucleons to larger momenta up to order {1.5-2)ppn at sat- 
uration density. Descriptions of higher density neutron 
matter with softer interactions if they are consistently 
evolved to lower scales, must include 3n (and potentially 
4n) interactions. 

Phenomenological two nucleon potentials such as the 
Argonne potential have been constructed to describe 
scattering data up to relative momenta ~ 600 MeV with 
high accuracy [13j . Despite the fact that the Argonne po- 
tential has been fit up to laboratory energies of 350 MeV, 
it very well reproduces scattering data up to much larger 
energies [T3] The AV8' interaction we employ in this 
study is identical to the full AV18 interaction in s and p 
waves, and includes the dominant onc-pion interaction in 
higher partial waves. Chiral interactions also reproduce 
the scattering data very well below 350 MeV laboratory 
energy, but they fail rapidly above because of the cutoff 
in presently available interactions. At larger momentum 
transfer, the potentials cannot describe inelasticities, but 
in scattering channels where inelasticities are known to be 
small they have been shown to provide a good descrip- 
tion. They also provide good predictions [IS] of high- 
momentum components of nuclear wave functions as ob- 
served in nucleon [inilTT] and electron scattering [181ll9j . 
These high momentum observables provide a test of the 
assumed short-distance features. In the low-energy high- 
momentum region relevant to neutron stars the inelastic- 
ities in 2n scattering must be absorbed into many-body 
forces (3n, 4n, . . .) intimately connected to the short- 
distance behavior of the 2n interaction. 

The nuclear Hamiltonians we consider contain the non 
relativistic kinetic energy, and the 2n and 3n interactions: 

H = -^+V2n+V3n- (1) 

For the 2n potential, we use the Argonne AV8' model [20] 
and the form of the 3n interaction is inspired by both 
the Urbana IX and the Illinois models [3]. We consider 
a range of 3n interactions that contain long-distance s 
and p wave 27r exchange contributions, an intermediate- 
range (37r loops) contribution, and a spin-independent 
short-range repulsive term. Explicitly, 

(2) 

This form of interaction includes all the terms present in 
low order chiral interaction, plus selected terms found to 
be important in studies of light nuclei and nuclear matter 
using the Argonne interactions. 



The structure of the operators O appearing above are 
defined in Ref. [5]. The relative contributions of these 
four components of the 3n force depends on the 2n inter- 
action. We find that for the Argonne potential, the 2n 
interactions suppress the long-distance (27r) contribution 
of the 3n force in the ground state. This suppression is 
a result of the pion-range correlations induced by the 2n 
force, we find it also occurs for the super-soft core NN in- 
teraction [5T] . For typical ranges of values of the strength 
parameters and considered in Ref. [5] we find 

the contribution of these operators to the ground state 
energy is repulsive but very small at all densities studied. 
In contrast, this interaction is large and attractive in light 
nuclei where both neutrons and protons contribute. The 
intermediate-range (37r) 3n interaction was introduced to 
fit the properties of weakly bound neutron-rich nuclei 
such as ^He 5J. Earlier calculations [10] have shown that 
this interaction is strong and attractive in neutron matter 
for typical values of A^-^ quoted in Ref. ^ . In this work, 
we explored a range of values for A^^^ from zero to that 
in the Illinois-7 3n interaction [22] because the structure 
of this term is still not fully understood or constrained. 
We use a phenomenological short-range repulsive term 
as in the Urbana and Illinois three-body forces, with 
Vr = ArO"' = ARj2cycT^im.n,)T^m^rjk), where 
the function T{x) is defined in Ref. |5j. We have also 
considered a different form = ^Rj2cyc'^i^ij)'^i'''jk) 
with and v{r) — exp(— 2/^r); other different forms of Vr 
have been explored, giving very similar results. 

The 3n interaction we employ is not intended to be a 
microscopic treatment of the complete 3n interaction. It 
assumes that for the neutron matter equation of state 
the effects of more complicated spin-dependent short- 
distance 3n interactions, relativistic effects, and potential 
4n interactions can be mimicked with simplified three- 
neutron interactions with a wide range of spatial depen- 
dence. This assumption has been tested in the case of 
relativistic corrections, where in Ref. [23j it was found 
that the density dependence of the relativistic effects is 
similar to that of the 3n interaction. Further tests of the 
density dependence of specific higher-order terms in the 
chiral interaction are valuable. The different forms of Vr 
we have explored span a wide range of density depen- 
dence for the 3n interaction, as shown below. 

For the 3n interaction we vary both A^^,- and p to study 
the sensitivity to short-range physics. The strength of 
the short-range 3n interaction Ar is taken to be a free 
parameter adjusted to yield the experimentally acces- 
sible nuclear symmetry energy. Although not proven, 
we make the following reasonable assumptions: (1) rel- 
ativistic effects in neutron matter show a similar den- 
sity dependence to the short-range three-nucleon inter- 
action as carefully studied in Ref. [53], (2) the density 
dependence of additional spin-dependent short-range 3n 
interactions (for example, higher-order terms in chiral ex- 
pansions) in the equation of state of neutron matter can 
be described in a spin-independent model, and (3) four- 
nucleon force contributions with different density depen- 
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FIG. 1. (Color online) The energy per particle of neutron 
matter for different values of the nuclear symmetry energy 
(-Esym). For each value of i?sym the corresponding band shows 
the effect of different spatial and spin structures of the three- 
neutron interaction. The inset shows the linear correlation 
between -Bsym and its density derivative L. 



dence are suppressed relative to the 3n force for densities 
up to (2 — 3)/9o- This last assumption can be justified at 
nuclear density by the high-precision fits to light-nuclei 
obtained with only 3n forces [53] ; at higher density this 
model assumption can be tested by its predicted correla- 
tion between properties of neutron-rich nuclei and neu- 
tron stars. 

We assume that E^sym = -Bnoutron(/5o) - i^nuclcar(/5o) 

and using experimental values of Esym = 32 ± 2 MeV 
pS] and -Enucioar(Po) = —16.0 ± 0.1 MeV from nuclear 
masses models we obtain an empirical constraint for 
neutron matter energy -Encution(po) = 16 ± 2 MeV. Po- 
tential higher-order corrections to the quadratic nuclear 
symmetry energy, for which there is some theoretical mo- 
tivation but no clear experimental evidence, may affect 
the extraction of the neutron matter energy and increase 
the associated error. In this work we ignore these poorly 
known corrections and tune to reproduce the neu- 
tron matter energy in the range 16 ± 2 MeV. Our results 
are shown in Fig. [T] where the green and blue points 
are QMC results for different choices of A/j correspond- 
ing to £:„outro„(Po) = 16 MeV ( E^^^ = 32 MeV) and 
£^ncutro„(po) = 17.7 McV ( -Bsym = 33.7 MeV), respec- 
tively. The results are compared to those obtained using 
a 2n force without 3n (i?sym = 30.5 MeV), and 2n com- 
bined with the Urbana IX 3n {E^y^a = 35.1 MeV). The 
bands depict the sensitivity to short-distance spin and 
spatial structure of the 3n interaction and are obtained 
by varying the range of the 3n short-distance force and 

In the vicinity of nuclear density, i?ncutron(p) = 
£'ncutron(po) + L/'i [p — po) / Po wherc L is related to the 
derivative of the nuclear symmetry energy. The inset in 
Fig. [T] shows the correlation between E'sym and L. This 
correlation is insensitive to the large variations in the 



range of the short-range 3n force fi and the strength 
of the 37r term ^3^. This is in sharp contrast to the 
predictions of mean field theories where the slope was 
found to be very sensitive to the choice of effective in- 
teractions [27 \ . Previous calculations of neutron matter 
up to pn [28j use a chiral 2n interaction fit to laboratory 
energies of 350 MeV plus the two-pion exchange three- 
nucleon interaction to calculate the neutron matter equa- 
tion of state using perturbation theory. In contrast to 
our results, a significant repulsion from the 27r exchange 
long-range 3n interaction was found. Since this force is 
better constrained by light nuclei, these earlier calcula- 
tions can make a prediction for the neutron matter energy 
independent of the phenomenological short-range inter- 
action, which plays an important role in our calculation. 
To understand this basic difference, further tests of the 
convergence of perturbation theory and the chiral expan- 
sion in the diagrammatic calculations, a survey of other 
two-body interactions in the AFDMC, and the incorpo- 
ration of chiral interactions in non-perturbative methods 
such as lattice and suitable extension of QMC would be 
necessary. 

Current determinations of L have relied on analysis 
of neutron-skins, surface contributions to the symme- 
try energy of neutron-rich nuclei, and isospin diffusion 
in heavy-ion reactions. These studies have been useful, 
but not very constraining as acceptable values are in the 
range L = 40 - 100 MeV [5S]. However, a better deter- 
mination of L even with modest reduction in the error 
would test our model for 2n and 3n interactions. 

The predictions of QMC can be accurately fit using 

E{p)=a (^J^y + b (3) 

where the coefficients a and a are sensitive to the low 
density behavior of the EoS, while b and /3 are sensitive 
to the high density physics [S^. We find that the 3n 
force plays a key role in determining the coefficient b and 
the variation of the other EoS parameters is compara- 
tively small. Numerical values for these parameters are 
reported in Table |l] for selected Hamiltonians. 
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2.49 


UIX 




35.1 


63.6 


13.4 


0.514 


5.62 


2.436 



TABLE I. Fitting parameters for the neutron matter EoS de- 
fined in Eq. [Sjfor selected different Ifamiltonians. 



To calculate the mass and radius of neutron stars we 
solve the Tolman-Oppenheimer-Volkoff (TOV) equations 
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FIG. 2. (Color online) Mass-radius relation for the EoS with 
three-neutron interactions corresponding to the bands for dif- 
ferent Ssym shown in Fig.jl] The intersections with the orange 
lines roughly indicate central densities realized in these stars. 



for the hydrostatic structure of a spherical non rotating 
star using the QMC equation of state for neutron matter 
[301 131] ■ The QMC EoS we use is for p > pcrust = 0.08 
fm"^. Below this density we use the EoS of the crust 
obtained in earlier works in Refs. and [33^. 

The neutron star mass-radius predictions are obtained 
by varying the 3n force and are shown in Fig. [2j The 
striking feature is the estimated error in the neutron star 
radius with a canonical mass of 1.4 Mgoiar- The uncer- 
tainty in the measured symmetry energy of ±2 MeV leads 
to an uncertainty of about 3 km for the radius, while the 
uncertainties in the short-distance structure of the 3n 
force predicts a radius uncertainty of ^ 1 km. The dif- 
ferent bands of Fig. [2] correspond to the EoS of Fig. [l] 
with the same colors, giving different values of i^sym- 
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FIG. 3. (Color online) Bounds on the maximum mass and 
radius for different equations of state as a function of the 
critical density pc. The left panel shows the maximum mass; 
the right top and bottom panels show the maximum possible 
radius for any neutron star with mass greater than 1.2A/soiar 
and for a neutron star with M — 1.4Msoiar, respectively. 

The central density of stars with A/^ l.SMgoiar are 
larger than 3po. At these higher densities, effects such as 
relativistic corrections to the kinetic energy, retardation 



in the potential, and four- and higher body forces become 
important. Consequently, non-relativistic models violate 
causality and predict a sound speed Cg = y^dp/de > c for 
p ~ (4 — 5)po- To overcome this deficiency we adopt the 
strategy suggested in Ref. [34j and replace the EoS above 
a critical density pc by the maximally stiff or causal EoS 
given by p{e) = c^e — ec, where p is the pressure, e is 
the energy density, c is the speed of light and Ec is a 
constant. This EoS is maximally stiff and predicts the 
most rapid increase of pressure with energy density with- 
out violating causality. The constant Ec is the parameter 
that determines the discontinuity in energy density be- 
tween the low- and high-density equations of state. Our 
choice of Ec ensures that the energy density is continuous 
and provides an upper bound on both the radius and the 
maximum mass of the neutron star. 

Figure [3] shows how the bounds on the maximum ra- 
dius and mass of the neutron star vary with our choice of 
the critical density pc- It also illustrates that the bounds 
provide useful constraints only when the EoS is known up 
to (2 — 3)po- In Ref. [35 bounds on the radius were de- 
rived by using an EoS of neutron matter calculated up to 
po with specific assumptions about polytropic equations 
of state at higher densities. Our upper bounds are model 
independent and show that the radius of a I.4Msoiar neu- 
tron star can be as large as 16 km if pc ^ po- To obtain a 
tighter bound the equation of state between Ipo and 2pQ 
is important. The red, green, blue and black curves are 
predictions corresponding to the 3n interaction strength 



fit to i?s- 



30.5, 32.0, 33.7 and 35.1 MeV, respectively. 



We also note that these bounds do not change much for 
Pc ^ 4po because the QMC EoS is already close to being 
maximally stiff in this region. These upper bounds pro- 
vide a direct relation between the experimentally measur- 
able nuclear symmetry energy and the maximum possible 
mass and radius of neutron stars. 

To summarize, we predict that the correlation between 
the symmetry energy and its derivative at nuclear den- 
sity is nearly independent of the detailed short-range 3n 
force once its strength is tuned to give a particular value 
of Esym- Consequently, in our model one short-distance 
parameter An completely determines the behavior of the 
EoS. At higher density, the sensitivity to short-distance 
behavior of the 3n interaction translates to an uncer- 
tainty of about 1 km for the neutron star radius with 
mass M — 1.4Msoiar- The uncertainty at high density 
due to a poorly constrained symmetry energy is larger, 
~ 3 km. Within our model we predict that neutron star 
radii are in the 10 — 13 km range for nuclear symmetry 
energy in the range 32 — 34 MeV. If nuclear experiments 
can determine that i?sym < 32 MeV, QMC predicts that 
L < 45 MeV at nuclear density, and for neutron stars it 
predicts M^ax < 2.2Msoiar and i? < 12 km for a neutron 
star with M — 1.4Msoiar ■ The relationship between the 
symmetry energy and its density dependence is exper- 
imentally relevant, and its implications on the neutron 
star mass radius relationship are subject to clear obser- 
vational tests. 
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